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Abstract: In this article we describe simulations of ZZ, W^Z and W~^W~ production 
based on the positive weight next-to-leading-order (NLO) matching scheme, POWHEG, in 
the HERWIG++ event generator. Building on earlier efforts within the HERWIG++ frame- 
work, the simulation includes a full description of truncated showering effects, required to 
correctly model soft, wide angle, emissions in angular-ordered parton showers. We utilize 
simple relations among each of the diboson cross sections, holding to 0{as), in constructing 
the simulation. Spin correlation effects are also included in the decays of the vector bosons 
at the tree order. A large part of this work is concerned with a full and thorough validation 
of the simulations through comparisons with alternative methods and calculations. 



Keywords: QCD, Monte Carlo, NLO Computations, Resummation, Collider Physics.. 



Contents 



|l]. Introduction Q 

H]. Hardest emission cross sections m 
|2.1| Kinematics and phase space ^ 



2.2 Differential cross section 



2.3 Matrix elements 10 



|3|. Event generation 12 

|3.1| Hardest radiation generation 12 

|3.2| Spin correlations and vector boson decays 13 

|3.3| Truncated and vetoed parton showers 13 



g. Results [14 
[4.1| Inclusive observables 15 



4.2 Exclusive observables 19 
|5]. Conclusion 26 
|A|. Regularized and unregularized splitting functions 28 



1. Introduction 

The production of electroweak gauge boson pairs is a subject of significant interest in 
collider physics. Weak boson pair production was first studied in detail at LEP2, at centre- 
of-mass energies of up to 209 GeV [1-5], with the aim of precisely measuring the W boson 
mass and with a view to directly probing and testing the non-Abelian character of the 
Standard Model gauge structure. In the period 1996-2000 approximately 40000 W^W~ 
and 1000 ZZ events were recorded by the four LEP experiments [6]. At the Tevatron 
W~^W~ and W^Z production have also been the focus of extensive analysis by the CDF 
and DO collaborations [7-15], albeit with an estimated factor of twenty less diboson events 
per experiment, based on their most recent publications. Recently the observation of ZZ 
production has also been reported by CDF and DO. 

Although all of these LEP2 and Tevatron measurements are in agreement with Standard 
Model predictions, it is generally considered that should non-Standard Model physics play a 
role in this process, its effects will be manifest greater at greater energy scales. In particular, 
the effects of anomalous WWj and WWZ triple gauge boson couplings, are predicted to 
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grow with the invariant mass of the gauge boson pair [16]. Thus, even though fewer events 
win be recorded, the high energy data stih being accumulated by CDF and DO combined, 
has the potential to surpass the sensitivity of the LEP2 data to such effects. At the LHC the 
inclusive cross section for vector boson pair production is around a factor of ten higher than 
at the Tevatron, moreover, the luminosity is also substantially higher. It is then expected 
that in the near future very large quantities of high energy weak boson pairs will be recorded 
by the LHC experiments and that the non-Abelian nature of the weak bosons will be tested 
with a new level of vigor. To facilitate such precision tests theoretical predictions, ideally in 
the form of fully exclusive Monte Carlo simulations, should be made as accurate as possible. 

Aside from being an interesting process in its own right, weak boson pair production 
deserves to be thoroughly understood given its role as an important background in many 
searches for currently undiscovered particles, most notably the Higgs boson. Notwithstand- 
ing the fact that the latest Tevatron analysis excludes, at the 95% confidence level, the 
Standard Model Higgs boson from having a mass in the range 158 < mn < 175 GeV [17], 
if the mass is above ~ 135 GeV the Higgs boson will primarily decay into or ZZ 

pairs [18-22]. Moreover, alternative models with extended Higgs sectors can yield large 
branching fractions for charged Higgs bosons decaying into W^Z pairs [23]. W^Z pair 
production is also a problematic background in studies of strong WW scattering [24] and 
in many supcrsjmnnetry search channels, specifically, those associated with a trilepton final 
state and missing transverse energy [25-27]. As with the dedicated studies of weak boson 
pair production discussed above, the discovery potential of these ongoing and imminent 
analyses depends heavily on the quality of the theoretical inputs. 

In general, the main source of uncertainty associated with leading order theoretical pre- 
dictions for hadronic collisions comes in the form of next-to-leading order QCD corrections. 
Next-to-leading order (NLO) QCD corrections to ZZ, W^Z and W^W~ production were 
first calculated and studied in the early nineties by two separate groups [28-33] . Some years 
later these results were improved upon by Dixon et al, who performed the calculations at 
the level of helicity amplitudes, thus providing full knowledge of the O (as) corrections to 
these processes including the leptonic decays of the massive vector bosons [34,35]. Fur- 
ther improvements were made by Campbell and Ellis who extended the results of the latter 
work beyond the narrow width approximation, including contributions from singly resonant 
Feynman diagrams and interference effects between intermediate Z bosons and photons [36] . 

Following these results, in the early part of the last decade, a number of ground breaking 
developments took place in the field of Monte Carlo event generator research, most signif- 
icantly, the invention of the CKKW(-L) and MLM algorithms, combining parton shower 
simulations together with those based on tree level matrix elements [37-42] and, separately, 
the Mc@NLO [43] and Powheg [44, 45] formalisms for consistently including fully differ- 
ential NLO corrections in parton shower simulations. In fact the Mc@NLO formalism was 
first successfully demonstrated through the simulation of I^+W^~ pair production and the 
first public release of the program comprised of W~^W^, W^Z and ZZ production, two 
years after the work of Campbell and Ellis. An impressive number of important Standard 
Model processes can now be simulated with the MCQNLO program [46-49], in addition, the 
modeling of unstable particle production has been enhanced through the inclusion of full 
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spin correlations in the leading order and real contributions by the method developed in 
Ref. [50]. Nowadays a similar number of processes may also be simulated with publicly 
available programs based on the alternative POWHEG formalism [51-59], first realized in 
the case of ZZ hadroproduction in 2006 [60]. 

Lately significant steps have been taken in the direction of automating the Mc@NLO 
and PoWHEG methods [59,61,62]. The most recent of these efforts features PoWHEG 
simulations of W^W~ and ZZ production, integrating the calculation of Campbell and 
Ellis within the framework of the Sherpa event generator [62]. This simulation there- 
fore includes singly resonant contributions, in particular Z/^* interference effects in W^Z 
and Z pair hadroproduction, which are not included in the program discussed here. For 
reasons of technical simplicity, greatly increased computational efficiency, and minded to 
best utilize the existing HERWIG++ infrastructure, e.g. the facility to include higher order 
QCD and QED corrections to the decays of vector bosons, we have based our work on the 
original calculations of Frixione et at [28,32,33], as in Mc@NLO, valid in the double pole 
approximation. 

We also note that many important collider physics analyses, such as the study of the di- 
vector boson signal [8,10,11,63-67] and anomalous triple gauge boson couplings [10,68,69], 
employ invariant mass cuts on the vector boson decay products - almost always to Z 
bosons - which generally reduces the singly resonant contributions to negligible levels. 
The same is true of Higgs boson searches involving final-states comprising of Z bosons 
[70]. In the case of the W boson the same invariant mass cuts around the boson decay 
products are not necessary since, in the case of the Z, the singly resonant contributions 
are only able to become sizable through Z/7 interference. Thus, in W pair production 
the singly resonant terms are subject to a simple Breit-Wigner suppression. Accordingly, 
the program presented here can be applied to the simulation of WW pair production for 
study as a signal process and as a background to Higgs boson production, as an alternative 
to the Mc@NLO program [13,71-73]. Furthermore, as is usually the case in higher order 
calculations involving unstable particles, an ad hoc scheme is adopted in Ref. [36] to include 
the width of the vector bosons which, in this case, mistreats the singly resonant parts, 
particularly in the doubly resonant region. Nevertheless, it is incumbent upon us to point 
out that should analyses not include some form of cut on Z boson decay products, limiting, 
in particular, Z/7 interference contributions, the matrix elements of Campbell and Ellis will 
continue to provide a good description of the production and decay processes, whereas those 
used here, like those of the underlying LO HERWIG++ simulation, will fail. Analyses not 
employing such cuts are typical of general searches for physics beyond the Standard Model 
e.g. the search for SUSY in the trilepton channel, for which singly resonant contributions 
are expected to add around 15% to the W^Z doubly resonant background ^. 

In what follows we present an application of the PoWHEG method to Monte Carlo 
simulations of weak boson pair production within the HERWIG++ event generator [74,75], 
including a full validation. In Sect. ^ we briefly recall the basic PoWHEG algorithm and 
give details concerning how the NLO cross section is organized in a form amenable to 

study of the impact of singly resonant diagrams on such an analysis can be found in Ref. [36]. 
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its implementation; the section closes with a discussion of the associated NLO matrix 
elements [28,32,33], noting exact relations between the W^Z, W~^W~ , and ZZ next-to- 
leading order cross sections. In Section ^ we describe the key points in the simulation 
process, including the implementation of truncated and vetoed parton showers, occurring 
after the hardest emission has been generated. In Section ^ we show results from our 
implementation, comparing it to two independent NLO programs, MCFM and MCQNLO, 
before summarizing our findings in Sect. ^. 

2. Hardest emission cross sections 

The starting point for all POWHEG simulations is the so-called hardest emission cross 
section, specifically, a matching between the constituents of the exact NLO differential 
cross section and the corresponding leading-log resummed cross section implicit in parton 
shower simulations [44]. For simple processes, such as the one we are considering, it can be 
written simply as 

da = B d^B 

where $b are Born variables, which fully determine the kinematics of leading order con- 
figurations, and <I>ij are radiative variables, parametrizing the kinematics of the hardest 
emission with respect to -B and B {^b) are the leading and next-to- leading order 
differential cross sections respectively, while R{^b,^r) is simply the bare, tree level, real 
emission cross section. The B function may be expressed as 

B = B {^b) + V + I d^RR {^B, ^r) , (2.2) 

V {^b) being the finite combination of, unresolvable, soft emission and virtual loop correc- 
tions, while R{^B, ^r) corresponds to the remaining, regularized, real emission corrections. 
The PoWHEG Sudakov form factor for the hardest emission, (pt), is defined as 

(pt) = exp 

where {^b, ^r) tends to the transverse momentum of the emitted parton in the soft and 
collinear limits. Emissions for which Jct < kT,min are considered as being unresolvable. In 
the ensuing subsections we give details concerning the components of the hardest emission 
cross section for the case of weak boson pair production (Eg. p.l|) , focusing on the B (^b) 
function and its regularization. 

2.1 Kinematics and phase space 

In this subsection we begin by giving details regarding our parametrization of the kinematics 
for the hadroproduction of a pair of weak vector bosons, with and without the emission 
of an additional radiated parton. For both classes of event we denote the momenta of the 
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partons incident from the ±2; directions by p@, while those of the weak bosons are labelled 
pi and p2- In discussing three-body, real emission, contributions the momentum of the 
additional, radiated, parton is denoted k. Also, since we shall frequently refer to the sum 
of the weak boson momenta, we define p = pi + p2. 

The parametrization of our two- and three-body kinematics is taken to be identical to 
that in Refs. [28,32,33]. To this end we first define precisely what we mean by the rest 
frame of the vector boson pair. In the context of three-body events, for qq and qg collisions, 
we shall use this term to refer to the frame in which the vector bosons are balanced in 
their three-momenta, with the incoming quark defining the +z axis and the transverse 
momentum of the other initial-state parton defining the +y axis; for gq collisions the gluon 
replaces the quark in defining the +z axis. For genuine two-body events the latter criterion, 
concerning the definition of the +y axis, is omitted since for such events and Pq are 
naturally acolincar. 

We now introduce a set of Born variables $_b = {p^, y, ^} and a set of radiative 
variables = {x, y, ip}, clearly defined as follows: 

p^ - the invariant mass squared of the vector boson pair 

y - the rapidity of the weak boson pair in the lab frame 

9 - the polar angle of pi in the rest frame of the vector boson pair 

(f) - the azimuthal angle of pi in the rest frame of the vector boson pair 

X - the ratio p'^/s where, as usual, s = (p^ + Pe)'^ 

y - the cosine of the polar angle of momentum k in the partonic 



Both and (f> range between and tt. Given a set of Born variables one can readily construct 
the momenta in the lab frame for two-body pg)+PQ — )• pi +P2 reactions, while augmenting 
these with a set of radiative variables enables one to fully reconstruct radiative processes 
p,^ + Pq Pi + P2 + k in the lab frame. For the latter case explicit expressions for the 
momenta in the rest frame of the vector boson pair are given in Ref. [32]. Since these are 
unwieldy, yet straightforward, we do not reproduce them here. When reconstructing the 
three-body events using the aforesaid formulae, all particles are then returned to the lab 
frame by applying the following Lorentz transformation, T, to each one: 



The first component of this transformation, M, is a rotation of angle arctan pt/\/p^ in the 

y — z plane, where pT is the transverse momentum of the vector boson system in the lab 
and partonic centre-of-mass frames: 



centre-of-mass frame 



T = ]B||]B_lM. 



(2.4) 
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Following R a transverse Lorentz boost, Bj^, is carried out such that the momentum of 
the stationary vector boson system, p, becomes {Et,0, —pT,0), where we have defined 



Et = yp'^+p'^- Lastly a longitudinal boost, By, gives the vector boson system rapidity 
y, returning all particles to the lab frame and completing the momentum reconstruction. 
Note that the ultimate step in the generation of all PoWHEG events involves randomizing 
the azimuthal orientation. 

The final kinematic quantities we wish to declare are the momentum fractions of the 
incoming partons with respect to the beam particles. These are not independent degrees 
of freedom but functions of #b and, in the case of three-body final-states, ^r. Taking P@ 
to represent the momenta of the parent beam hadrons, the momentum fractions x@ are 
defined according to the relation p@ = x^P@, whereupon it follows that 





^® = ^\l n .^l-x){l±y) ^"""^ ^@ = \/-e ^ (2-6) 

with s being the liadronic centre-of-mass energy. Note that for genuine two-body final 
states, namely those corresponding to leading order and virtual contributions, the limit 
X — >■ 1 is clearly implied in the evaluation of the all kinematics quantities, including the 
momentum fractions i. e. for such final states x@ = x@ . 

Having described the parameterization of the kinematics we move to specify the inte- 
gration measures in the two- and three-body phase spaces, for non-radiative and radiative 
events respectively. With these definitions in hand the phase space for the leading-order 
process can then be written as 

d$B = dx(sdxQd^B = ^dp^dyd^B, (2.7) 

where d^B the two-body phase space measure for the vector boson system. In the centre-of- 
mass frame of p, using conventional dimensional regularization, with n = 4 — 2e dimensions 

where p* is the magnitude of the three momentum of either weak boson in their rest frame. 

To parametrize the three-body phase space we factorize it into a product of two two- 
body phase spaces for the reactions p^ + Pq —>^p + k and p Pi + P2, by inserting the 
identity in the form of an integral over a delta function defining p as pi + p2- In this way 
we may readily write the phase space measure as 

d$ = d$B d^R / cr J {x, y) J (</.) , (2.9) 

(47r) x'^ \P J 

where^ 

J (x, y) = 22^ x^ (1 - x)'-^' (1 - y') , J (0) = lllyl sin-2^ cj> , (2.10) 



^An irrelevant overall factor ^ r(i-e) been dropped in writing J" {x, y) since it is equal to 1 + O (e^) . 
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and 



d$R = — dydxd0. (2.11) 
zvr 

The constant cr appears due to the use of dimensional regularization, it is given by 



cr = (47r) 



r(l-e)^r(l + e) 
r(l-26) 



(2.12) 



We emphasize that both 6 and (j) range between and vr, hence J d^R = 1. It is also worth 
noting that ^ f d(f> J {(j)) = I and lime^o J {4>) = ^■ 

Since we restrict ourselves to processes for which the NLO corrections contain at most 
soft and initial-state collinear singularities, the product of with the squared real emission 
matrix elements will be finite throughout the radiative phase space. With this in mind, 
following Refs. [28,76], we extract a factor of p"^ from J' (x, y) and then expand it in powers 
of e, to give 



J{x,y) = [SS{l-x) + C{x) (25{l + y) + 26{l-y))+'H{x,y)] ^ 

s 



in which 
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(2.14) 
(2.15) 

(2.16) 



where 7] = \f\^^ and p = (ra\ + 771-2)^ /s, ra\ and being the masses of the weak vector 
bosons. The /o-distributions appearing in Eqs. 2.15 and 2.16 are defined according to the 
relation 



dx h (x) 



In" (1 - x) 
1 — X 



1 In^fl-x) 

dx (/i(x) -/i(l)) ^ ^, (2.17) 

1 — X 



for any sufficiently regular test function h (x): in this case, the product of with the real 
emission matrix elements. 

For completeness, we note that the kinematic boundaries are 



(mi + 771,2) < < s , 
X (y) < X < 1 , 



|y|<-^in(^ 
|y| < 1, 



(2.18) 



with X (y) given by 
X (y) = max 



1+y^l-y 

X© O Xq 



(2.19) 



-7- 



In order to ease the numerical implementation of the p distributions we map the x variable 
into X, defined according to 

x{x,y) = x{y) + fi{yf X , f)[y) = -x{y). (2.20) 

Whereas the x integration domain was dependent on y, the x integral simply ranges from 
to 1. 

Finally we wish to clarify that the expressions for J {x,y) in Eqs. 2.1C1| and 2.13 are 
equivalent only up to terms of O (e); these terms do not contribute to the differential cross 
section in the limit e — t- 0. 

2.2 Differential cross section 

In this section we enumerate the various contributions to the NLO differential cross section. 
These are obtained by simply considering the product of the matrix elements with the phase 
space measures as written in Eq. |2.7| and Eqs. 2^^.17| , exploiting simplifications arising in 



the soft (j; — )• 1) and collinear limits (y — )• ±1) to integrate out radiative variables. The 
discussion here is rather similar to that in our earlier work [56], including some minor 
changes and clarifications, indeed the formulae in that article were derived with the current 
application in mind, so we shall be brief. 

In the following we shall refer to the parton types incident from the + and —z directions 
as a and b respectively. We denote the parton distribution function (PDF) for a parton 
of type i inside a beam liadron traveling in the ±z direction by f® {x(^,fi'jp), where fip is 
the factorization scale. For brevity we then introduce the luminosity function. Cab, as the 
product of the PDFs associated to a and b: 

Cab (^e, ^e) = /a® {x(B,f^F) (^e, /«f) • (2.21) 

The leading order contribution to the differential cross section is given by the product 
of the leading order spin and colour averaged squared amplitude, A^^, together with the 
luminosity function and fiux factor: 

da''' = B{<PB)d^B, B{'^B) = ^MZ{<^B)Cab{x^,xe) . (2.22) 

In all cases, due to the universal nature of the infrared divergences in the one-loop ampli- 
tudes, the virtual corrections to Mab generically written in the form 

M^l {<^b) = Vo M^, {^b) + Mir i^B) , (2.23) 
where the first term is universal and divergent, with Vq given by 



r 2 2 tt^ 



' g2 ^P"- "-9 



(2.24) 



2tt Vp^ 

while the second is process dependent and finite as e — )■ 0. For completeness, we define 
the bare virtual cross section and Vq {^b) function by analogy to the corresponding leading 
order quantities: 

1 

2^2 



= ^0 , ^0 ($b) = ^ Mil i^B) Cab (xe, xq) . (2.25) 
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The differential cross sections for the real emission processes, a + b 
following general form 



n + c, take the 



A^afe i'^B, '^r) Cab {X(B,Xe) d$ . 



2s 



(2.26) 



For each flavour combination we consider these corrections to comprise of three compo- 
nents, corresponding to the three terms 5, C (x) and Ti {x, y) in the phase space Jacobian, 
Eqs. 2.14- 2.1^ . We shall refer to these components as the so/i, collinear, and hard / resolved 
contributions to the cross section. 

The squared matrix elements for the real emission processes, in which a gluon is emitted 
from an initial-state quark or antiquark, factorize in the limit that the gluon is soft (x — )• 1) 
according to 

lim A^ffe i^B, ^r) = SirasfJ.''' -\2CfM^^, . (2.27) 

Thus, for X = 1 the integrand in Eq. p.26 is entirely independent of the radiative phase 
space and it becomes trivial to integrate over <l>/j. In doing this one finds the following 
expression for the soft contribution to the differential cross section: 



ascr ( n 



2-K 



Cf 



In7?-M61n% B{<^b) d<^B- 



(2.28) 



2 vr" 
? ~ T e 

All other sources of real corrections involve the emission of a quark or antiquark from 
external initial-state gluons, as such they do not contribute to the cross section in the soft 
limit; by contrast to the case of gluon emission, the matrix elements for such processes are 
regular in the limit x — ?■ 1, hence when they multiply the factor of p'^ in J' (x, y) (Eq. |2.13" ) 
the result is zero. 

We shall now turn our attention to the collinear limits. Since the real emission cor- 
rections to the process we are considering do not involve internal gluon lines, in the limit 
±1 the spin averaged squared matrix elements factorize trivially according to^ 

1 



y 



J/-S>±1 



['^B, ^r) = SVTQS^ 



2e 



(l-x)P„:-(x;6) M!^,{^b) 



xp\ 



(2.29) 



where i = a or b for y = ±1 respectively and c denotes the type of the emitted parton. 
Considering only those terms in the phase space Jacobian proportional to C (x) we then 
obtain the following expression for the collinear contribution to the cross section 



Co® 



da' 



(x) 

tic 



c® 



27r 

as \ 
2tt X 



p- 



a-c (P.ic + 41nr,) Q + In (^^)) B (ci>^) d^^, 



C®. (x) C% {X, 
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/i^X 
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Xq) B{^b) d^Bdx 
In (1 — x) 



(2.30) 



(1 



1 



X 



(l-x)P.^(x)-P.'^ (x), 



B (<I>b) d<l>Bdx, 



'in Eqs. 2.27 and 2.29| we show only the leading term in since all others ultimately vanish due to the 



factor of in Eq. 2.13. 
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where i = a in the case that parton a sphts to produce parton c, and i = b for the 
case that parton b branches to produce c. The functions P.t-^{x;€) are the bare Altarehi- 
Parisi sphtting kernels in n = 4 — 2e dimensions, exphcit expressions for which can be 
found throughout the hterature e.g. Refs. [56,77], while P.~ (x) denotes their regularized 
counterparts (Appendix^); the soft-collinear cross section, da^^®, is entirely due to the 
6 {I — x) term present in the latter. Lastly we declare £^ to be the ratio of the luminosity 
function evaluated at y = ±1 with respect to that found in the leading order cross section 
viz 



J/=±l Lab [X^, Xq) 



^Tb (^ffi' ~ '^af' (^ffi' , , ' ^ab (^0) ^g) " ~ — ' (2.31) 



The singular parts of the collinear contribution proportional to the regularized Altarelli- 
Parisi functions, da^^®, are exactly canceled by collinear counterterms in the MS scheme, 
thus they play no further role in our discussion. 

Finally, considering the third part of the J {x,y) phase space Jacobian, we have that 
the hard / resolved contribution to the real emission cross section is given by 



d^ab = ? - Cab (^e, Xq) B b) R , (2.32) 

ZTT X 



^ab = x,l-^ '^lf;^^\^> n{x,y). (2.33) 



Combining all of the various pieces together we are able to write the next-to-leading 
order differential cross section as 

da = B ($b) d^B + V ($b) d^B + R (^B, ^R) d^B d^R, (2.34) 

where the second term is given by the finite sum of the bare virtual corrections, d(T^°, 
the soft real emission contributions, do"^", and the soft-collinear contributions, do"^'^®. 
The i?($B,$ij) function is comprised of the remaining collinear and hard / resolved real 
emission corrections, 

R ($B, $^) = - IZab Cab {x^.xq) B ($s) , (2.35) 
zvr X 

IZab = 2C®~ {x) 5{l-y) + 2Cf~ (x) 5{l + y)+ Uab ■ 

Although it is not explicitly stated above, a summation over all contributing channels and 
flavour combinations is implied. Furthermore, in the case of the qg and gq initiated real 
emission corrections, one of the terms proportional to C®~ ix) in the TZab function should 

iic 

be omitted, either 2C®~, (x) 6 {1 — y) if 6 is a gluon, or 2C®~ (x) 6 {1 + y) if a is a gluon, 
since collinear branchings of the initial-state quarks do not occur in those processes. 

2.3 Matrix elements 

As noted in the introduction the simulation described in this article uses the NLO QCD 
calculations for ZZ, W^Z and W~^W~ production of Frixione et al. [28,32,33]. We have 
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observed that these computations share a common underlying structure which we have 
exploited here. In particular we find that all W^Z production matrix elements, leading- 
order, real and virtual, are simply related to the corresponding ZZ production matrix 
elements by the following replacements: 



ez , mw mz , 

where the quantities on the left are written using the notation of the W^Z production 
publication [32] and those on the right use that of the ZZ article [28]. In the latter work 
g'd,L and g^^i^ denote the left handed couplings of the Z boson to up- and down-type quarks, 
Fij is the W boson coupling to quark flavours i and j, while ez is the trilinear gauge coupling 
(ez = 9u,L — gd,h)- In Ref. [28] and are the vector and axial- vector couplings of the 
Z bosons to the colliding quarks. 

We also flnd that the W^Z matrix elements are simply related to those of W~^W~ 
production by 

K^.^l, el ^^(p2-my^f (p2), 

m,L ^ , - "^w) 4' {P') , (2.37) 

9ifl. ^ mz ^ mw , 



where on the left Kij denotes the relevant Cabibbo-Kobayashi-Maskawa matrix element in 
W^Z production, while on the right i is used to refer to the flavour of the colliding quarks 
(up- or down-type). The functions c|* (p^) and c** (p^) are the coefficients of those parts 
of the W~^W~ matrix elements corresponding to s-channel trilinear gauge-boson graphs 
interfering with themselves and with t-channel graphs respectively. 

The validity of these relations was examined analytically using Mathematical. In our 
simulation we have employed the matrix elements of Ref. [32] for W^Z production and 



applied the transformations in Eqs. 2.36 and 2.37 to these when generating ZZ and W^W 



production events respectively. The correctness of relations Eqs. 2.36 and 2.37 at NLO is 
tested again, numerically, by comparing to alternative calculations in MCFM and MCONLO. 

We attribute the fact that these relations hold at the NLO level to the deceptively simple 
Dirac structure in the real and virtual corrections. In the beginnings of Refs. [28,32,33] it is 
observed that, for both sets of radiative corrections, all of the Dirac traces can be expressed 
in the form Tr [(a + 675) F], where a and b are constants and F is an arbitrary string of 7 
matrices, excluding 75. In the case of the virtual corrections the 75 term never contributes 
to the final trace: it leads to a term proportional to the antisymmetric Levi-Civita tensor, 
for which there are not enough linearly independent momenta available to give rise to a 



*A Mathematica file detailing these checks for all cross section formulae in Refs. [28,32,33] is available 
from the author on request. These checks also reveal that, in both cases, all matrix elements satisfy 
additional symmetries involving the momenta of the final state vector bosons before the last transformation, 
mw — >■ "T-z / mz — >■ mw, is carried out; further exploration of this point is beyond the scope of this work. 
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non-vanishing contribution. Similarly, in the case of the real corrections this term gives a 
purely imaginary contribution, where a real part must be taken. In fact, it is remarked in 
Refs. [32, 33] that the NLO W^Z and W~^W~ cross sections could be considered to have 
arisen from a fictions theory, with suitably redefined vector couplings and no axial ones. 

Although the results in [28,32,33] were updated in Ref. [34,35], at the end of the event 
generation process we simulate the decays of the vector bosons according to either the full 
2 — )• 4 or 2 — >■ 5 tree order matrix elements, using basically the same procedure as is adopted 
in several other Mc@NLO and POWHEG simulations [47,48,50,51,57,78]. The 2—^4 and 
2 — )• 5 particle helicity amplitudes were constructed using the C++ libraries present in 
ThePEG event generator tool kit [79], based on the Helas formalism [80]; in practice we 
have implemented these matrix elements as a convolution of the 2—7-2 and 2—7-3 particle 
production spin density matrices with 1—7-2 particle decay matrices (Sect. |3.2D . We have 
checked that the corresponding 2 — )• 2 and 2 — )• 3 'undecayed' matrix elements reproduce 
the results obtained with those in Refs. [28,32,33], however, since the latter are considerably 
faster to evaluate we refrain from using them until the last step of the hardest emission 
event generation i.e. the vector boson decays. Further details of the decay simulation 
procedure are deferred to Section ^. 



3. Event generation 

In this section we describe the steps by which the event simulation is carried out in practice: 
the generation of the hardest emission about a 2 — )- 2 underlying Born configuration, the 
decays of the vector boson pairs according to 2 — t- 5 particle matrix elements, and the 
subsequent parton showering. 



3.1 Hardest radiation generation 

The first step in the simulation process involves generating a set of Born variables, ^b, and 



hence a 2 — )- 2 kinematic configuration, according to the B {^b) function (Eqs. |2.2| and |2. 341) . 



The B (^b) function is numerically implemented directly according to the formulae given 



in Section g, having applied the transformation in Eq. 2.20 to x. To this end we generate a 
set of Born and radiative variables {^b,^r} by sampling Eq. 2.34 using a VEGAS based 
algorithm [79], the ^b. a-re then simply discarded, leaving ^b distributed according to the 
integral with respect to ^b, in other words B (^b)- 

With $5 in hand we proceed to generate the hardest radiation according to the square 



bracketed term in Eq. ^T|. The exponent in the PoWHEG Sudakov form factor Eq. 2.3 
consists of an integral over a sum of terms, 

B \^b) 27r X 

one for each real emission process, ^Lab being equal to Hah with the plus and p regularization 
prescriptions omitted. Here, in implementing the generation of we have opted not 
to generate x and y directly but rather we re-express the d$R integration measure in 
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terms of pT and y^, where is rapidity of the hardest emission in the hadronic centre-of- 
mass system. Making this change of variables removes the sHghtly awkward ^-function in 



Eq. 2.3, replacing it by a lower bound on the integration over px- The distribution of the 



radiative variables viz. the square bracketed term in Eq. 2.1, is then sampled using the veto 
algorithm [81]. 

With the resulting set of Born and radiative variables {^b,^r} we may then fully 
reconstruct the kinematics of the 2 — )• 3 hardest emission events directly using the general 
formulae given in Ref. [32]. Alternatively in the rare event that pT < kT,min, where in 
this work we have chosen kT^min = 2 GeV, the emission is considered to be unresolvable, in 
which case only the 2 — )• 2 kinematics are reconstructed before proceeding to the next step 
of the event generation. 

Finally we note that in generating for the hardest emission we have used px as 
the factorization scale in the PDFs and the renormalization scale in the strong coupling 
constant, in accordance with the DDT formulation of the Sudakov form factor [60, 82]. 
This completes the generation of the hardest emission kinematics according to Equation 



2.1 



3.2 Spin correlations and vector boson decays 

Having generated a set of 2 — )• 3 kinematics we now calculate and store the production 
spin density matrix, ■^^^\^X2\2^ from the associated tree order helicity amplitudes, where 
the pairs of indices {Ai,A2} and {Ai,A2} label the helicities of each vector boson in the 
amplitudes and conjugate amplitudes respectively. These helicity amplitudes were con- 
structed using the C++ libraries present in ThePEG event generator tool kit [79], based 
on the Helas formalism [80]. A set of two-body decay kinematics are then generated for 
each vector boson, isotropic in their rest frames, from which corresponding decay matrices 
M':\ andA^r are calculated and contracted with the production matrix element. The 

AiAi A2A2 

decay kinematics are then kept provided 



where 7?. is a random number in the range [0, 1], while and M^"^ correspond to the 

traces of the decay matrices. If the decay kinematics are rejected the process is repeated, 
using newly generated sets of momenta for the decay products, until the inequality in 



Eq. 3.2 is satisfied. By generating the decays of the vector bosons in this way, breaking the 
2 — )• 5 process down into a 2 — )• 3 process followed by two-body decays, we avoid the more 
intensive operation of computing the full 2 — ?■ 5 body helicity amplitudes and summing over 
helicity amplitudes repeatedly. Lastly, we note that for the tiny fraction of events in which 
Pt < kT,min the Same procedure is applied to generate the decay kinematics, albeit using 
the helicity amplitudes of the 2 — )• 2 leading order process to compute the production spin 
density matrix as opposed to those of the 2 — )• 3 processes. 

3.3 Truncated and vetoed parton showers 

In order to shower the hardest emission configurations, we first compute a set of parton 
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shower branching kinematics, <J>^ , corresponding to the hard radiation. More specif- 
ically, we precisely determine the HERWIG++ branching variables [83] which would have 
been assigned to exactly this momentum configuration had it been generated by initiating 
the parton shower from the underlying Born configuration^. Having ascertained how the 
hardest emission event may be reproduced by the usual parton shower apparatus, we return 
to the underlying Born configuration and proceed as follows: 

• the external leg deemed to have produced the hardest emission is evolved from the 
default shower starting scale to that of <I>^^^^, with the imposition that intervening 
branchings conserve flavour and have transverse momenta less than px- the truncated 
shower [44]. 

• the set of branching parameters ^^'^'^ is then inserted in this shower. 

• the evolution continues down to the cut-off scale, vetoing any emissions with trans- 
verse momenta greater than pT- the vetoed shower. 

• the non-emitting leg is evolved from the default shower starting scale down to the 
cut-off scale with a further vetoed shower. 

In the event that the hardest emission occurs in a region of phase space inaccessible to the 

parton shower, i.e. the wide angle / high pT dead zone [83], subsequent emissions will have 
sufficient resolving power to see the widely separated emitters individually. It follows that 
no truncated shower is then required, since this models coherent, large angle emission from 
more collimated configurations of partons, and so we proceed directly to the vetoed shower. 



4. Results 

In the following we present predictions from our POWHEG simulations in comparison with 
results obtained by independent calculations and alternative methods. The main aim of this 
work is to provide a robust validation of the simulations, such that they may be considered 
suitable for use in real physics analyses. 

All of the attendant tests have been carried out at nominal Tevatron and LHC centre- 
of-mass energies, respectively, 1.96 and 14 TeV. In all programs we have elected to use the 
MRST2002 NLO PDF set [85] interfaced through the LHAPDF package [86]. To analyze 
jet structure in the events we have used the /c_L-jet measure with the R parameter set 
to 0.7 [87,88], as implemented in the FastJet jet finder package [89], to carry out the 
associated clustering. 

^For technical details and formulae pertaining to this inverse reshuffling procedure we refer the reader 
to Refs. [54,84]. 
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4.1 Inclusive observables 



In order to check the calculation of the POWHEG differential cross section and B (^b) func- 
tions, Eqs. 2.1-2.2, we have compared our predictions for total cross sections and numerous 
inclusive observables against those of the NLO Monte Carlo calculator MCFM [36,90,91]. 
Since MCFM computes these quantities in fixed order perturbation theory and since we 
wish to test the various components of the simulation systematically, the PoWHEG results 
shown here have been obtained at the parton level, prior to showering. In other words, 
the predictions from our simulations here solely reflect the implementation of the hardest 
emission cross section, Eq. 2.1. 



To facilitate these investigations we have chosen to work with a fixed value of 100 GeV 
for the renormalization and factorization scales in both MCFM and the B{^b) functions 
in the PoWHEG simulations. Again, in order to have a meaningful comparison with MCFM 
we have run it in a mode where the narrow width approximation is assumed, moreover, we 
have arranged for each program to use the same vector boson decay channels. 

On the left of Fig. |l|we show the invariant mass of the weak boson system, p^, one of the 
Born variables for this process, for which there is excellent agreement between the PoWHEG 
result and that of MCFM. Predictions for the y Born variable, on the right-hand side of 
Fig. 1^, are also seemingly indistinguishable from the corresponding MCFM results. Recall 
that in the PoWHEG framework the Born variables are exactly preserved in the process 
of generating the hardest emission: they are distributed purely according to the B{^b) 
function. Consequently, if the PoWHEG simulations and the underlying NLO calculations 
are implemented correctly these Born variables must agree exactly with fixed order NLO 
predictions. The results in Figure |l] already therefore constitute a very sensitive test of our 
implement at ion . 

The distribution of the third Born variable 9 is strongly reflected in Figure ^, which 
shows the polar angle between the incident parton traveling in the +z direction and the 
first of the produced vector bosons in their rest frame; in W^Z and W^W^ production 
these are taken to be the and bosons respectively. This quantity is different to 
the 6 Born variable in that the latter is defined as the polar angle of the first vector boson 
with respect to the quark in qq and qg collisions, and the gluon in gq collisions, which may 
or may not be traveling in the +z direction in the lab. Nevertheless, since this variable 
is fully inclusive and a close relative of the Born variable, 9, the level of agreement shown 
here between the PoWHEG and MCFM predictions provides further strong assurance as to 
the correctness of our implementation. We add that we have also examined rapidity and 
pseudorapidity distributions of the individual vector bosons from our programs (not shown) 
which, like those in Figures || and ^, exhibit no discernible deviations from those given by 
MCFM. 
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Figure 1: In this figure we show predictions for the invariant mass (p^) and rapidity (y) of the 
vector boson pair system, in the left- and right-hand columns respectively; the results obtained 
using the Powheg simulation are shown in red while the blue dotted line and the black points 
represent the leading and next-to-leading order predictions from Mcfm. Since and y are Born 
variables in the Powheg simulation, they are distributed purely according to the B function. 



hence they must follow exactly the corresponding NLO prediction (Sect. 2.1). 
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Figure 2: The polar angle between the incident parton traveling in the +z direction and the 
first of the produced vector bosons in their rest frame; in W^Z and W^W~ production these are 
taken to be the and bosons respectively. Since this variable is fully inclusive and a close 
relative of the Born variable, 9, the level of agreement shown here between the Powheg and Mcfm 
predictions provides strong confirmation as to the correctness of our implementation. 



In Figure |3| we have displayed a number of distributions sensitive to the details of 
the decays of the vector bosons, specifically, the polar angle of one of the leptons pro- 
duced by one of the decaying, resonant, vector bosons in its rest frame and, separately, 
the corresponding transverse momentum spectra. In all cases the agreement between our 
predictions and those of MCFM is remarkably good. It is interesting to note that although 
our simulation only includes spin correlation effects in the leading order and real emission 
contributions to the cross section, it nevertheless reproduces very well the distributions pre- 
dicted by Mcfm, which also includes the effects of NLO virtual corrections at the level of 
the spin correlations. Corroborating evidence for such behaviour can be found in the work 
of Grazzini [92,93], whose calculations dealt with full NLO corrections to spin correlation 
matrices together with soft gluon resummation effects. 
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Figure 3: On the left we show the distribution of the polar angle of one of the leptons emitted by 
one of the decaying weak bosons, in the rest frame of the decay, while on the right hand-side we 
show the corresponding transverse momentum spectrum. The colouring of the different predictions 
is as in the previous figures. Note that in plotting these quantities the branching fractions of the 
vector boson decays have been divided out. 
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Actually, at Tevatron energies we find that the leading order production spin density 
matrix, calculated using the underlying Born kinematics {^b), yields more-or-less identical 
distributions to those shown in Figure ^ on the contrary, at the LHC this somewhat naive 
procedure produces cos 6i distributions close in shape to the leading order prediction, which 
is markedly different to the NLO one. These observations are very much in keeping with 
others in the literature [50], suggesting that virtual corrections to spin correlation matrices 
are typically small, whereas real emission corrections (which are generally larger at the LHC 
than at the Tevatron) can lead to sizeable effects. 

Finally we wish to draw attention to the tails of the lepton pT distributions. In Fig- 
ure |3| we see that for distributions forecast at Tevatron energies the POWHEG and MCFM 
predictions agree very well over the whole spectrum, whereas at the LHC we see that the 
two distributions overlap identically in the low px region but in the high pj- region the 
PoWHEG result is approximately 30% above that of MCFM. 

In fact this level of disagreement with respect to fixed order NLO predictions, in this 
region of phase space, is not unexpected. Events in which the final-state leptons have very 
high transverse momenta will, by their nature, typically contain associated high energy 
QCD radiation; this fact is substantiated by the shape of the corresponding leading order 
predictions shown in the blue dotted lines. Recall that, in the PoWHEG hardest emission 
cross section, the term responsible for the generation of the radiative variables is not simply 
equal to the NLO real emission cross section, but rather the NLO real emission cross section 
multiplied by a factor B (^b) /B {^b) and the Sudakov form factor. In the high pT limit 
the latter factor tends to one, hence, PoWHEG events with high transverse momentum are 
generated according to the NLO real emission cross section multiplied by B (^b) / B {^b)- 
Loosely speaking this factor is characteristic of the NLO total cross section if -factor, thus 
one expects the rate of events including high pT emissions to be different in PoWHEG 
with respect to a pure NLO calculation by such an amount (different by terms beyond 
NLO accuracy). This is indeed what we observe in Figure g Finally we add that the 
Born variables are unique exceptions to this reasoning since they are fully preserved, by 
construction, when radiation is generated, hence, unlike other inclusive observables the Born 
variables will always exactly equal the NLO prediction, regardless of whether the events are 
associated to high px emissions or not. 

4.2 Exclusive observables 

In this section we shift the focus of our validation onto observables which more directly 
assess the generation of the hardest emission (Sect. p.lD and additional radiation arising 
from the subsequent vetoed and truncated showers (Sect. |3.3| ). To this end we compare 
our results to two different approaches, namely, the default angular-ordered parton shower 
simulation in HERWIG++ and also Mc@NLO v3.4 [94]. Since the latter program matches 
the same NLO matrix elements with the older Herwig parton shower [95,96] it is formally 
of equivalent accuracy to our PoWHEG simulation, on the other hand, the former includes 
only LO matrix elements, hence, it is anticipated that it will fail to adequately model high 
Pt radiation. 
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Figure 4: The transverse momentum spectrum of the produced weak boson pair system at Teva- 
tron (left) and LHC (right) energies. Predictions from the Mc@NLO and HERWIG++ PoWHEG 
simulations are present as black and red dashed lines respectively. Results from the leading order 
HerwigH — h parton shower simulation are also shown as blue dotted lines. For the case of a single 
emission this quantity is equivalent to the radiative variable px introduced in Section 2.1. 



All of the results presented in this subsection, from each of the three approaches, were 
obtained at the parton level, after parton showering. The predictions from our PoWHEG 
simulation within HERWIG++ are displayed as red dashed lines, while those of Mc@NLO 
and HERWIG++ with the PoWHEG feature disabled are shown as black and blue dotted 
lines respectively. 

Figures ^ and ^ exhibit the transverse momentum spectra of the di- vector boson system 
and the hardest emitted jet. These distributions directly reflect the px dependence of the 
hardest emission cross section Eq. 2A, modulo small smearing effects due to the truncated 
and vetoed parton showers. All three approaches are seen to agree well in the low px regions, 
where the parton shower approximation is expected to be reliable. In the high px regions 
one can see that in all cases the parton shower by itself underestimates the production rate 
with respect to MCQNLO and PoWHEG, moreover, we note that the degree to which it is 
underestimated is worse at LHC energies, than at the Tevatron, where the available phase 
space is more constrained (c/. Fig. |^). 

In general the degree of overlap between the Mc@NLO and PoWHEG results is quite 
good, with differences between the two sets of predictions only being noticeable in the 
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Figure 5: The transverse momentum of the hardest jet in weak boson pair production at the LHC 
(left) and Tevatron (right), assuming a nominal LHC centre-of-mass energy, ^/s — 14TeV. The 
colouring of the histograms is the same as in figure ^. As for the case of weak boson pair system, 
for a single emission this quantity is equal to the radiative variable px in Section 2.1. 



high pt tails for LHC energies. Such behavior has already been observed, and its nature 
well documented, in publications concerning other POWHEG simulations [55,56]; it is ba- 
sically a further manifestation of the effect discussed at the end of Section [4.1| , namely, 
that the distribution of the hardest emission in the PoWHEG method (Eq. |2.lD is given by 
B ($_b) /B ($b) multiplied by the NLO real emission cross section, whereas in MCQNLO it 
is given by the real emission cross section alone. It follows that the PoWHEG predictions 
tend to exceed those of Mc@NLO and fixed order calculations when px is large, by a factor 
of the order of the NLO total cross section iT-factor. In keeping with this one sees that the 
relative differences seen in the tails of the transverse momentum spectra reflect the size of 
the relevant iC-factors (c/. Fig. ||) and, accordingly, they are somewhat larger at the LHC 
than at the Tevatron. 

Note that within the PoWHEG formalism it is possible to introduce so-called damp- 
ing factors [55], which act in such a way as to reduce the effects of the multiplicative 
B {^b) /B {^b) term, leading to cross sections in the high px domain closer to those of 
fixed order calculations. However, there is no theoretical motivation to implement, or not 
to implement such damping, since the contribution from this factor, in the regions of phase 
space corresponding to high pT emission, is formally subleading - the differences with re- 
spect to Mc@NLO and fixed order predictions being indicative of the associated theoretical 
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Figure 6: The transverse momentum spectra of individual weak gauge bosons in di-vector boson 
production at the LHC, assuming a hadronic centre-of-mass energy of 14 TeV. Predictions from the 
McQnlo and HerwigH — h POWHEG simulations are present as black and red dashed lines respec- 
tively. Predictions from the default HerwigH — h simulation, with the Powheg feature disabled, 
are represented by blue dotted lines. 



uncertainties. 

In Figure ^ we plot the transverse momentum spectra of individual vector bosons in 
each of the weak boson pair production channels, at nominal LHC energies only. In contrast 
to the transverse momentum spectrum of the diboson system, in fixed order perturbation 
theory this observable receives contributions at leading order for all values of pT- Thus 
one expects that any differences between the leading and next-to-leading parton shower 
predictions should be small, of order as- It is then remarkable that these distributions 
show the next-to- leading order parton shower predictions, in good agreement with one 
another, exceeding the leading order parton shower prediction by up to a factor of four at 
high PT- 

This curious result has already been noted and investigated in the original calculation 
of the fixed order NLO corrections to W^Z production in Ref. [32]. The same effect 
was subsequently observed in the case of W'^W~ production and subject to the same 
analysis in Ref. [33]. The detailed studies carried out in these publications conclude that 
the enhancement seen at high pT is greatly dominated by contributions arising from qg 
initiated real emission corrections. The reason for this large qg contribution was considered 
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to be twofold. Firstly, it was noted that the luminosity for the qg channel was more 
than one order of magnitude greater than that of the qq channel at the LHC; secondly, 
in qg reactions, when the radiated quark and one of the weak bosons are produced with 
sufficiently high transverse momenta, the other weak boson may be produced as a 'soft' 
emission from the recoiling quark - a process which carries a large logarithmic enhancement 
factor, log^ {pt z/''^\v^ ■ further noted in Ref. [33] that fewer partonic subprocesses 

can participate in this enhancement mechanism in the case of W^W~ production than in 
W^Z production, accordingly, we observe that the magnitude of the effect is somewhat 
smaller in the latter case. 

Having studied several pj- spectra we now turn to examine other distributions sensitive 
to the generation of the radiative variables and subsequent parton showering. As noted in 
earlier works concerning Higgs boson production, Ref. [56], the rapidity correlation between 
the leading jet and the recoiling colourless system is an interesting observable to examine 
from this point of view: for the case of a single hard emission the rapidity correlation yfe — y 
can be expressed purely in terms of In order to provide some additional physical insight 
regarding the nature of this quantity, we note that in the limit where the radiated parton 
is produced in the region perpendicular to the colliding beam partons, in the partonic 
centre-of-mass frame 

where 6^ denotes the polar angle of the emitted parton in that frame. Furthermore, when 
the radiated parton is emitted along the ibz directions, yfc — y tends to ±00. 

In Figure ^ we show predictions for yu — Y distributions in ZZ production and Z 
production at the Tevatron and LHC respectively. For each process we have considered how 
the results are affected by varying the pT cut on the leading jet. The general trends seen 
in these plots are qualitatively the same as those obtained in Higgs boson production and 
Higgs boson production in association with a W^/Z boson [56], so too are our conclusions 
relating to them. 

We remind the reader that, in general, parton shower Monte Carlo programs may 
not populate the full real emission phase space, resulting in so-called dead-zones. This is 
certainly true of the Herwig and HERWIG++ simulations. The presence of dead zones 
in the real emission phase space follows directly from the scale choice used to initiate the 
parton shower evolution^, moreover, they are typically located in regions of phase space 
associated with high px emissions. In Herwig and HERWIG++ the dead zone for the 
first emission in processes such as this, comprised of a single initial-state colour dipole, is 
centred on 0^, = ^ (yk— y = 0), moreover, the angular breadth of the dead zone increases 
symmetrically and monotonically about this point with the energy of the emitted radiation 
(see e.g. Fig. 7 of Ref. [56]). The dip feature seen in the ordinary parton shower results 
(blue dots) directly reflects the angular characteristics of this unpopulated region of phase 
space: in all cases, as the px cut on the leading jet increases the dip becomes broader and 
deeper. 

®For explicit phase space computations and maps concerning the origin of dead zones and their connection 
to the choice of the initial evolution scales see Refs. [56,83,97]. 
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Figure 7: In this figure we plot tlie difference in rapidity between tlie liardest jet and tire di- vector 
boson system, with different cuts imposed on the transverse momentum of the leading jet. On the 
left hand side we show predictions for these observables in ZZ production at the Tevatron, while 
on the right hand side they are shown for the case of Z production at the LHC. 



Since Mc@NLO and PoWHEG aim to fully include NLO corrections within the parton 
shower simulation, they naturally populate all of the real emission phase space. In keeping 
with this, we observe that the predictions from these two approaches do not show the same 
significant dip in the central region of the distributions. Whereas the PoWHEG simulation 
fills this phase space independently of the detailed workings of the parton shower to which it 
is subsequently attached, the Mc@NLO approach involves carefully augmenting the parton 
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shower simulation by the difference between its own approximate real emission cross section 
and the true real emission cross section in the NLO calculation. In particular, this means 
that the distribution of radiation from McONLO in the dead zone follows exactly the fixed 
order NLO calculation, while either side of it the distribution differs from this by O (ag) 
terms. 

Having noted this NNLO discontinuity in the radiation pattern, it is then understand- 
able that the MCONLO predictions (black) can exhibit some minor irregularities and dif- 
ferences with respect to POWHEG (red) in the central region and that these should reflect, 
somewhat, the trends seen in the results obtained using the parton shower alone. On the 
contrary, the response of the PoWHEG predictions to the increasing pj- cut on the leading 
jet lends itself to a more straightforward interpretation based on simple kinematic reason- 
ing, namely, that the y^ — y distribution should become more central, as the phase space 
available for small angle emissions - which populate the tails - becomes reduced relative 
to that available for large angle emissions. Furthermore, we point out that for observables 
employing cuts which exclude the softer regions of phase space, such as y^ — y with a px 
cut of 80 GeV on the leading jet, one expects that the PoWHEG predictions exceed those of 
conventional NLO calculations due to the B {^b) /B ($b) factor multiplying the real emis- 
sion cross section in Eq. 2A . This behaviour is apparent in Figure 0. Finally, we remark 
that we have reproduced all of these distributions using our PoWHEG simulation with the 
truncated shower feature disabled, with no observable consequences. 

Figure |8| shows the differences in rapidity between the leading jet and one of the vector 
bosons in W~^W~ production at the Tevatron and W''Z production at the LHC. Since 
these distributions are closely related to those in Fig. ^, we argue that the differences seen 
in the leading order parton shower predictions with respect to MCQNLO and PoWHEG are 
again attributable to the dead zone in the former; a fortiori considering the variation of the 
uncorrected parton shower results with respect to the changing pT cut on the leading jet. 
We contend that the peculiar shape of the pure parton shower results reflect the impression 
left by the dead zone in y^ — y convoluted with the rapidity distribution of the vector bosons 
with respect to one another. In all of the distributions the agreement between the MCONLO 
and PoWHEG results is quite satisfactory: note that both approaches formally only offer a 
leading order description of this quantity. Some small distortion can be seen on the right 
of the Mc@NLO distribution at the Tevatron, for a pT cut of 80 GeV on the leading jet, 
which we tentatively suggest is indicative of the asymmetric parton shower prediction. 

LHC predictions for jet multiplicity distributions in W^W~ and W~ Z production can 
be found in Figure ^, assuming a hadronic centre-of-mass energy of 14 TeV. As one expects, 
in all cases jet multiplicities decrease rapidly as the pj- cut on the leading jet is increased. 
One can also see that the results obtained using the parton shower alone show a tendency 
to overestimate the number of events without any jets in comparison to the other methods. 
This is also expected given our earlier discussions concerning the dead zone in the parton 
shower phase space. Finally we note that the PoWHEG predictions for the number of events 
with one jet are above those of MCONLO by an amount characteristic of the NLO i^-factor. 
As noted previously, this systematic effect can be directly attributed to the presence of the 
B (^b) /B {^b) factor multiplying the real emission part of the PoWHEG hardest emission 
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Figure 8: On tlie left of this figure we sliow tlie difference in rapidity between the liardest jet and 
the W~ boson in W^W~ production at the Tevatron, with pr cuts of 10, 40 and 80 GeV imposed 
on the transverse momentum of the leading jet. Analogously, on the right hand side, we show the 
rapidity correlation between the hardest jet and the Z boson in W~ Z production events. 



cross section (Eg. Once again we note that the presence of this term modifies the 

distribution of hard radiation by terms of NNLO significance only. 



5. Conclusion 

In this article we have presented an implementation of the PoWHEG NLO matching formal- 
ism for simulations of weak boson pair production and decay, in the double pole approxi- 
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Figure 9: Here we sliow jet multiplicity distributions in W^W^ and W^Z production at the 
LHC; given a centre-of-mass energy of 14TeV. The three pairs of distributions, running from the 
top to the bottom of the figure, result from applying three different sets of cuts (10, 40 and 80 
GeV) to all of the jets in each event. 



mation. These simulations have been integrated within the HERWIG++ Monte Carlo event 
generator, including truncated shower effects to account for colour coherence phenomena. 

In constructing this NLO event generator we have employed novel relations between 
the W^Z cross sections and those of W^W~ and ZZ production. Total cross sections and 
parton level NLO distributions were found to be in excellent agreement with predictions 
obtained from the the MCFM NLO Monte Carlo program. 
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The shapes of the emission spectra from the full simulation, including parton shower 
effects, are seen to generally compare well with those of MCONLO in a wide variety of 
kinematic distributions - both of which exhibit large corrections with respect to the default 
parton shower predictions. Where minor differences have arisen between our results and 
MCONLO they have been studied in detail. 

As noted in previous works comparing Mc@NLO and POWHEG [55, 56] we observe 
a tendency for PoWHEG to produce slightly more hard radiation than Mc@NLO. The 
explanation given for this effect in those publications is seen to hold well here, specifically, 
that the B(^b) /B{^b) factor which multiplies the real part of the PoWHEG hardest 
emission cross section leads to an enhancement of high px radiation with respect to the 
corresponding NLO prediction; the differences being formally of order ag. As in Refs. [55,56] 
we also find that the Mc@NLO program exhibits a sensitivity to the phase space partitioning 
in the underlying parton shower simulation (c/. Figs. ^ and however, as with the 
enhancement of hard radiation in PoWHEG, this is formally representative of NNLO effects. 

All weak boson boson pair production simulations presented here are due for inclusion 
in the next public release of the HERWIG++ event generator. 
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A. Regularized and unregularized splitting functions 



We write the 'customary' regularized Altarelli-Parisi functions in terms of p-distributions 
as 

P.^ (x) = PP- {x)+C.^ (p.^ +4.lnr]] 5(1- x) , 

I.IC ^ ' i ir.^ ' I.IC \^I.IC I I ^ ' ^ 



where 
pPg{x) = 2CA 

PI;,{x) = Cf 



X 1 — X ,^ , 

H \- x( \ — x) 



{l-x) 
l + x^ 



l + il-x) 



PPg (X) = Tr \x^ + (1 - XY 

and 



Cgg — Ca, 
Cqq = Cf, 



Pgg 
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_ 3 

Pqq — 2' 



bn = — ( — Ca Tnnu 

47r V 3 3 



with all other p. ~ and C ~ being equal to zero. 
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